Automatic Filters for the Detection of Coherent Structure in Spatiotemporal Systems 
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Most current methods for identifying coherent structures in spatially-extended systems rely on 
prior information about the form which those structures take. Here we present two new approaches 
to automatically filter the changing configurations of spatial dynamical systems and extract coherent 
structures. One, local sensitivity filtering, is a modification of the local Lyapunov exponent approach 
suitable to cellular automata and other discrete spatial systems. The other, local statistical com- 
plexity filtering, calculates the amount of information needed for optimal prediction of the system's 
behavior in the vicinity of a given point. By examining the changing spatiotemporal distributions 
of these quantities, we can find the coherent structures in a variety of pattern-forming cellular au- 
tomata, without needing to guess or postulate the form of that structure. We apply both filters 
to elementary and cyclical cellular automata (EGA and CCA) and find that they readily identify 
particles, domains and other more complicated structures. We compare the results from EC A with 
earlier ones based upon the theory of formal languages, and the results from CCA with a more 
traditional approach based on an order parameter and free energy. While sensitivity and statistical 
complexity are equally adept at uncovering structure, they are based on different system properties 
(dynamical and probabilistic, respectively), and provide complementary information. 
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I. INTRODUCTION 

Coherent structures are ubiquitous in nature, the re- 
sult of complex patterns of interaction between simple 
units P, . Such structures are not simply epiphenom- 
ena of the microscopic physics, but instead often govern 
the system's macroscopic properties, providing powerful 
collective degrees of freedom — they make good "handles" 
on the system 3j. Condensed matter physics provides 
a host of equililDrium systems in which this is true Q, 
for instance in antiferromagnets 5] and liquid crystals 
[i^, and the anomalous properties of the high tempera- 
ture superconductors may be due, in part, to the pres- 
ence of ordered charge density waves (stripes) 7]. Co- 
herent structures also govern the behavior of many non- 
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equilibrium systems Q^. The obvious examples are bio- 
logical systems, which are rife with coherent structure, 
and are most certainly far from equilibrium (until they 
die) H, . Another example may be found in economics 
where the spatial pattern of economic activity has been 
postulated to be governed by the emergence of coher- 
ent structures in the form of cities and hierarchies of ur- 
ban centers 0. It has even been proposed that coher- 
ent structures can perform adaptive information process- 
ing tasks, and are therefore the typical physical embodi- 
ment of emergent computation (l^, llJ? HH • This 
conjecture is supported by both simulations of biologi- 
cal morphogenesis and experiments on the adaptive 
responses of stomata in plants 17]. In view of the in- 
trinsic scientific interest of complex, spatially-extended 
systems and the potential power of coherent structures 
to describe them concisely, it is important to develop 
techniques which can uncover and define coherent struc- 
tures in both the equilibrium and non-equilibrium cases. 
In this paper, we present two such filters, applicable to 
spatially extended discrete systems. 

Most existing methods for identifying coherent struc- 
tures rely on knowing some details about the structure 
one is looking for. Matched filters are designed so that 
signals of certain known character maximize the response 
[iq . Such filters may match either the structures them- 
selves (e.g., 19j), or the uninteresting details of the 
background in which the structures are embedded (e.g., 
poL [2l( ) . The latter technique has been particularly use- 
ful in equilibrium systems with broken symmetry. The 



2 



identification of an appropriate order parameter and an 
associated free energy implicitly defines a filter which 
matches the background, i.e., the ordered state. Depar- 
tures from the ordered state, such as domain walls, vor- 
tices and other topological defects constitute mesoscopic 
structure. While extremely successful, this method re- 
quires a lot of trial and error and some knowledge of the 
underlying micro-dynamics. In far-from-equilibrium sys- 
tems, much of the progress on detecting coherent struc- 
tures has involved the application of regular language 
theory [23 to one-dimensional deterministic cellular au- 
tomata [23, 13 1^ • Specific filters, adapted to par- 
ticular cellular automata rules, have been designed to de- 
tect structures dynamically jll HHHIIO, El • In these 
cases as well however, one must know a lot about the sys- 
tem's dynamics to work out, by hand, a suitable filter. 
Furthermore, to apply such notions to cellular automata 
of higher dimension, one would have to extend regular 
language theory to such dimensions, which presents sig- 
nificant difficulties . Extensions to stochastic systems 
would be even more problematic. 

The two filters introduced in this paper address these 
issues. Both are "automatic" filters, in that they require 
no prior information about the system's microscopic dy- 
namics or the structures generated. They can be applied 
to both equilibrium and non-equilibrium systems of, in 
principle, any dimensionality. One of them is applica- 
ble to systems with stochastic dynamics. The two filters 
are, however, based upon very different system proper- 
ties, m introduces the local sensitivity^ an adaptation 
of Lyapunov exponents to cellular automata, measuring 
the degree to which small perturbations can alter the dy- 
namics. ^Illl defines the local statistical complexity, which 
measures the spatiotemporal distribution of the amount 
of information needed to optimally predict the system's 
dynamics. It is therefore a probabilistic measure, com- 
pletely compatible with any underlying stochastic dy- 
namics. Both methods identify coherent structures by 
filtering the system's changing configuration with respect 
to sensitivity (or complexity), and tracking the chang- 
ing spatial distributions of these fields. We apply each 
filter to the task of detecting coherent structures in cel- 
lular automata. Sections ^11 Bl and ^III CI look at "el- 
ementary" (one-dimensional, binary, range-one) cellular 
automata (EGA), where we find that both filters read- 
ily identify EGA particles, domains and domain walls, 
although they emphasize these structures to differing de- 
grees. In section ^IVI we use cyclic cellular automata 
(GGA), a self-organizing model of excitable media, to 
compare the local sensitivity and statistical complexity 
with a more traditional order parameter approach. (We 
identify the appropriate order parameter for this system 
in ^IV Ah Our conclusion gives a summary and discusses 
some open questions. Finally, appendices discuss previ- 
ous attempts to adapt the Lyapunov exponent to cellular 
automata (0, and provide details on the construction of 
the GGA order parameter and effective free energy (|B)) . 



II. LOCAL SENSITIVITY 

The first of our filters for coherent structures is based 
upon a measure of local sensitivity: the instability of the 
discrete cellular automata field under local perturbation. 
Filtering with respect to local sensitivity is motivated 
by a desire to distinguish autonomous objects from the 
rest of the GA's configuration field. By "autonomous 
objects" we mean those structures that have the great- 
est influence upon the future dynamics of the system. 
Perturbing an autonomous object will affect the future 
a great deal. In contrast, perturbations of the rest of 
the system, the dependent parts, should be quickly over- 
ridden and "healed" by the autonomous dynamics. The 
set of autonomous objects constitutes a high-level model 
of the system, and knowledge of them should allow us 
to infer most of the rest of the dynamics, as well as the 
objects' own futures. 

A famous example of a system governed by au- 
tonomous objects is rule 110 of the elementary cellular 
automata. A rule 110 configuration field (Fig.Q) tends to 
exhibit particles, and knowledge of the position of these 
particles carries considerable predictive power [2^ [s^ . A 
second example may be found in the cyclic cellular au- 
tomata (see below, ^IV|) . For certain parameter values 
the GGA configuration field evolves into competing spi- 
ral waves, the cores of which are autonomous (see Figure 
0). 




FIG. 1: (Color online.) Typical space-time diagram of rule 
110; time advances from left to right. Note the presence of 
the regular background domain, and the particles which stand 
out in the background by contrast. 



For these two examples, rule 110 and spiral GGA, the 
presence of structure is obvious. However this is not al- 
ways the case (see, e.g., Ref. 30j, and the discussion of 
rule 146 in MBI below). Even when it is, reliably iden- 
tifying all the structure present in a GA's configuration 
field is difficult to do by eye and it is often not at all 
clear how autonomous such structures are. One might 
therefore hope for an automatic filter capable of not only 
distinguishing particles, domains, and extended objects 
like strings and domain walls in large spatiotemporal sys- 
tems, but also inferring their degree of autonomy. 



3 




FIG. 2: (Color online.) Phenomenology of the cyclic CA 
(defined in gTv]) in one of its spiral- forming phases (a^ = 4, T = 
2, r = 1). Time progresses from top to bottom, the panels 
depicting (a) the random initial conditions, (b) the formation 
of spiral wave cores, and (c) their eventual domination of the 
entire lattice. 

A. Defining Local Sensitivity 

Since we want to filter cellular automata with respect 
to autonomy, a natural point of departure is some math- 
ematical formalization of the idea that a system's future 
may be very sensitive to details of its present state. De- 



vaney's definition of chaos |3J, and the notions of to" 
logical and metric (Kolmogorov- Sinai) entropy rates j 




are two ways in which this intuition can be formalized. 
Here, however, we will begin with the idea of a dynam- 
ical system's Lyapunov exponents, since these provide 
a more refined and quantitative measure of sensitivity- 
dependence than either the mere fact of chaos, or than 
the metric entropy rate (which is generally the sum of 
the positive Lyapunov exponents). We briefly review the 
definition of Lyapunov exponents, before developing our 
filter. 

A one-dimensional map /, for which the system's state 
is given by a single real variable has a single Lyapunov 
exponent, given by the following limit: 

A= lim -log|D/"(xo)| 

where D is the derivative operator. (The limit can be 
shown to exist and be identical for almost all initial points 
xq.) The interpretation is that, if one begins with two 
points xq and xq + e, the distance between their future 
images, after n time steps, is approximately ee^^; this 
relation becomes exact as e ^ 0. In the case of a dynam- 
ical system described by m state variables, the Lyapunov 
exponents are defined as the m eigenvalues, Ai . . . Am, of 
the matrix 

A= lim (Dt/»(xo)Z)r(xo))'^'" 

n^co 

where D is again the derivative operator, and denotes 
the adjoint of A. These give the exponential growth rates 
of small perturbations applied along different directions 
(which are not, however, the eigenvectors of A, but those 
of Df). The spectrum of Lyapunov exponents thus gives 
a fairly detailed picture of the kinds and degrees of sen- 
sitivity displayed by the system dynamics. 

These definition take no account of the system's spa- 
tial structure; this should be rectified in an application to 
cellular automata. In addition, we must define the expo- 
nents in a way which is spatially locals so that their varia- 
tion over space can be studied, and structures identified. 
In continuous, spatially-extended systems the traditional 
approach has been to consider the global configurations 
of the system as points in a Hilbert space 35, p. 53], 
and then expand the above definition to include infinite- 
dimensional spaces. Although this produces a meaningful 
set of global exponents, these exponents contain no in- 
formation about whether the system is more sensitive at 
particular points in space. Such information is crucial for 
identifying coherent structure. Also crucial is an ability 
to identify structures at varying spatial scales. For ex- 
ample the elementary CA rule 110 has particles that can 
be over twenty cells wide — and all parts of the particle 
need not be equally sensitive. Some more recent work 
has attempted to define Lyapunov exponents in a spa- 
tially local manner, for precisely these reasons 0. 

A second challenge is to define Lyapunov exponents 
in a way which is applicable to systems which are dis- 
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cretized, both in state and space. The classical defini- 
tion of the Lyapunov exponent invokes limits taken as 
two initial conditions approach arbitrarily close to one 
another in continuous state spaces. (This is implicit in 
the use of derivatives.) In cellular automata, in which 
the state is discretized, the closest distinct configurations 
possible are two configurations differing only on one cell^ . 
Space is also discretized: there are no distances interme- 
diate between one and two differing cells. Several at- 
tempts have been made to frame definitions of Lyapunov 
exponents which accommodate these facts; we discuss 
these attempts, and why we feel they are not suitable 
for our present purposes in Appendix H Here we take 
the straightforward approach of perturbing small patches 
of contiguous cells, evolving the system forward in time, 
and measuring the distance between the perturbed and 
unperturbed configurations. This eliminates the need to 
take derivatives, and, by varying the size of the perturbed 
region, one can examine the structures present at differ- 
ent length scales. 

Formally, we define the local sensitivity at point r, to 
as the result of the following calculation. Let P be the 
set of all cells at distance at most p of r. The parameter 
p is called the perturbation range. Let S be the set of 
all possible configurations restricted to P (i.e., the set of 
words of length \P\ taking the states as alphabet). For 
one 5 G 5, we set the cells in P to their state in 5, and 
keep the original configuration for all other cells. Then 
we let the automaton evolve for / steps. We name the 
parameter / the future depth. Finally, we measure the 
area of the difference plumes, the total number of cells 
which differ from the original configuration, at each time 
step and compute the mean, weighting the area at time 
t by where w is the set of points at time t that 

depend on (r, to)- We average this on all s G 5 to get 
the sensitivity at that point, ^(r, to), which, intuitively, 
gauges how unstable the system is with respect to per- 
turbations there. Note that ^ is normalized so that it 
always lies between and 1, facilitating the comparison 
of distinct rules. 

Two parameters enter into our calculation of ^, though 
we suppress them in the notation: the perturbation range 
p and the future depth /. Both must be chosen with some 
care. If they are too small, the calculated ^ is excessively 
sensitive to fiuctuations inside the background domain or 
within particles, but over-large parameter choices tend 
to make things blurry, and are quite time-consuming. 
Choosing the right p is like adjusting a microscope to 
the size of the object one wants to observe. If the mag- 
nification factor is too high, one might miss the bigger 
structures, or even encounter diffraction artifacts (coun- 
terpart of being too close to the discretization scale). It 



There are metrics, such as the Cantor metric, in which CA con- 
figurations can approach arbitrarily close, but they are inappro- 
priate for local studies, such as ours. 



may seem self-defeating that we should need to know the 
size of the relevant structures to successfully use an auto- 
matic filter. However there may be structures at different 
scales, and varying p allows one to filter these preferen- 
tially. Similarly, the coherent structures of many systems 
can be arranged in a causal hierarchy, where higher-level 
objects drive lower-level ones, and this can be detected 
by increasing the future-depth /, which tends to pick 
out the most autonomous features. (See the example of 
EGA rule 146 below, especially Figure IHl obtained with 
/ = 30.) Note that the levels of the spatial and causal 
hierarchies need not be aligned with one another, i.e., the 
largest structures need not be the most autonomous. 



B. Results from Elementary Cellular Automata 

We now present the results of applying the dynami- 
cal sensitivity filter to several 1+lD elementary CA. In 
Figure 01 we show the configuration field of rule 110 (a) 
and the field filtered with respect to the sensitivity (6); 
Figure^does the same for rule 54. The filtering automat- 
ically distinguishes the autonomous features of the CA, 
in particular the particles' evolution through time. The 
autonomy of the particles is refiected by their dark tone. 
The domains, in contrast, are much less autonomous, and 
appear lighter. This makes intuitive sense. Domains are 
large, ordered regions of the CA, and their perturbation 
has little effect upon the dynamics, because any defects 
are either quickly healed or remain confined. Particles 
in contrast are relatively complex and localized objects 
which travel through the domains. Perturbing a parti- 
cle generally strongly changes the dynamics of the CA, 
either through its destruction, or via a significant alter- 
ation of its attributes — location, internal phase, type of 
particle, etc. 25]. 

The utility of the local sensitivity filter can be seen by 
comparing the raw configurations produced by rules 22 
and 146 to their filtered fields (Figures [51 and El respec- 
tively) . The configuration fields of the two rules look very 
similar, both appearing highly disordered. Once filtered, 
however, it is clear that they are quite different. Rule 22 
is chaotic, and all points have roughly equal, and strong, 
infiuence on the future of the system. All of this, along 
with the absence of autonomous objects, is plain from 
the sensitivity filter (Figure (S^). In contrast, rule 146 
(Figure O is composed of domains separated by bound- 
aries which wander over time. The domain walls appear 
as light traces, indicating that they are less autonomous 
than the domains they separate. The domain walls of rule 
146 are dependent objects because their motion is deter- 
mined by the chaotic dynamics of the domains to either 
side of them. This is a clear instance of the fact that 
while the presence of an autonomous object implies the 
existence of structures, the presence of structures does 
not necessarily imply autonomy. 

It is worth noting that in linear CA rules (those which 
are a sum (mod 2) on a subset of the neighborhood) 
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FIG. 3: Evolution of rule 110 from a random initial condi- 
tion (the same as that shown in Figure Q repeated for conve- 
nience): a, configuration field; b, sensitivity field (calculated 
with p = 1, / = 10); c, complexity field, using light-cones of 
depth 3. In b and c, higher values of the field are denoted 
by darker points. In c, the vertical grey band on the left is 
due to the fact that we need a past to compute the causal 
states. There are a lot of particles in the beginning, which 
makes the background domain rather unusual and complex, 
but it gradually bleaches. 



the local sensitivity is uniform over space. This is be- 
cause the Hamming distance between the original and 
perturbed configurations is independent of the perturbed 
cell, which in turn is because the difference plume always 
has the same shape. As expected from this argument, di- 
rect calculations of the sensitivity field on additive rules 
such as EGA 60, 85, 90, 105, 150 and 170 produce per- 
fectly uniform results. (We omit plots for reasons of vi- 
sual monotony.) 



III. LOCAL STATISTICAL COMPLEXITY 

Although the sensitivity is very effective at uncover- 
ing structures, as well as telling us how stable those 
structures are, it has a significant drawback. The GA 
configuration field must be actively perturbed many 




FIG. 4: Evolution of rule 54 from a random initial condition: 
a, configuration field; 6, local sensitivity (p = 1, / = 10); 
c, statistical complexity for rule 54 (light-cone depth of 3). 
In c, note that the background domain is light (low complex- 
ity) , the particle are grey, and most of the collisions are darker 
(highest complexity). The tiny particles that are merely phase 
shifts in the periodic background are made clear here, while 
they are hard to identify by eye in the configuration field 
(though misalignment between that field and your printer ma- 
trix can help!). 



times and the results compared to the unperturbed case. 
Such a procedure is computationally intensive and of- 
ten impossible to apply effectively to experimental data. 
We now discuss the local statistical complexity C{f^t) 
[IE 11^7 5 the calculation of which only requires obser- 
vations, rather than active and repeated perturbations. 
We shall see that the statistical complexity has other ad- 
vantages as well. 



A. Local Causal States and Their Complexity 

Let x(r, t) be an n + ID field, possibly stochastic, in 
which interactions between different space-time points 




FIG. 5: The evolution of rule 22 from a random initial con- 
dition: a, configuration field; 6, local sensitivity (calculated 
with p = 1, / = 10); c, statistical complexity. The nearly- 
uniform sensitivity field in b reflects the chaotic nature of the 
rule. 



propagate at speed c. As in ^39], define the past light 
cone of the space-time point (r, t) as all points which 
could influence x{f^t)^ i.e., all points {q^u) where u < t 
and 11^— r] I < c{t — u). The future light cone of (r, t) is 
the set of all points which could be influenced by what 
happens at (r, t). l~{f^ t) is the configuration of the field 
in the past light cone, and /^(r, t) the field in the future 
light cone. The distribution of future light cone configu- 
rations, given the configuration in the past, is P(/+|/~). 

Any function r] of l~ defines a local statistic. It sum- 
marizes the influence of all the space-time points which 
could affect what happens at (r, t). Such local statis- 
tics should tell us something about "what comes next," 
which is Information theory lets us quantify how in- 
formative different statistics are, and so guides our choice 
among them. 

The information about variable x in variable y is 



I[x;y] = (log2 



F{x)F{y) 



(1) 



where F{x^y) is joint probability, P(x) is marginal prob- 
ability, and (•) is expectation |4fl|. The information a 
statistic T] conveys about the future is /[/+; 77(/~)]. A 
statistic is sufficient if it is as informative as possible 
40], here if and only if /[/+;7/(r)] = /[/+;/"]. This is 
the same as requiring that P(/+|7^(/")) = P(/+|/"). 
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FIG. 6: The evolution of rule 146 from a random initial con- 
dition: a, the configuration field; 5, the sensitivity field (cal- 
culated using p = 1 and / = 30); c, the complexity field. 
Although the configuration field looks very similar to that 
of rule 22 (Figure |^ this field is not chaotic and has a do- 
main structure. Local sensitivity filtering (b) reveals the do- 
main walls as light (low sensitivity) traces bounded by darker 
zones. The domain walls aren't autonomous, their behavior 
is instead determined by what happens inside the domains, 
hence their low sensitivity. Note also the increase in sensitiv- 
ity when two walls are near: a small perturbation can lead 
them to merge (or prevent them from doing so). Statistical 
complexity filtering (c) reveals the domain walls as dark (high 
complexity) traces composed of localized triangular regions. 
Because the domain walls require more predictive informa- 
tion than the domains themselves, the statistical complexity 
is higher for the walls. The dark vertical lines are finite size 
effects and should be ignored 



A sufficient statistic retains all the predictive information 
in the data. Since we want optimal prediction, we confine 
ourselves to sufficient statistics. 

If we use a sufficient statistic r] for prediction, we must 
describe or encode it. Since r]{l~) is a function of /~, 
this encoding takes /[77(/~);/~] bits. If knowing rji lets 
us compute 7^2, which is also sufficient, then r]2 is a more 
concise summary, and /[7^i(/~); /~] > I[r]2{l~)]l~]. A 
minimal sufficient statistic AO] can be computed from 
any other sufficient statistic. In the present context, the 
minimal sufficient statistic is essentially unique (as dis- 
cussed below), and can be reached through the following 
construction. 

Take two past light cone configurations, l^ and 
Each has some conditional distribution over future light 
cone configurations, P(/+|/f ) and P(/+|/^) respectively. 
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The two past configurations are equivalent, ^ I2 ^ if 
those conditional distributions are equal. The set of con- 
figurations equivalent to /~ is [/"]. Our statistic is the 
function which maps past configurations to their equiva- 
lence classes: 

e{n^[l-] = {X : P(/+|A)=P(/+|r)} (2) 

Clearly, P(/+|e(r)) = P(/+|r), and so /[/+;e(/-)] = 
making e a sufficient statistic. The equiva- 
lence classes, the values e can take, are the causal states 
[371 141I \42L . Each causal state is a set of specific past 
light-cones, and all the cones it contains are equivalent, 
predicting the same possible futures with the same prob- 
abilities. Thus there is no advantage to subdividing the 
causal states, which are the coarsest set of predictively 
sufficient states. 

For any sufficient statistic r], P(/+|/~) = P(/+|7^(/~)). 
So if r]{l];) = r]{l^), then P(/+|/r) = PG+I^s"). and the 
two pasts belong to the same causal state. Since we can 
get the causal state from r]{l~), we can use the latter to 
compute e{l~). Thus, e is minimal. Moreover, e is the 
unique minimal sufficient statistic |3?, Corollary 3]: any 
other just relabels the same states. 

Because e is minimal, /[e(/~);/~] < I[r]{l~);l~], for 
any other sufficient statistic r]. Thus we can speak objec- 
tively about the minimal amount of information needed 
to predict the system, which is how much information 
about the past of the system is relevant to predicting its 
own dynamics. This quantity, /[e(/~); /~], is a character- 
istic of the system, and not of any particular model. We 
define the local statistical complexity as 



C{f,t) = I[e{l-{f,t));l-{f,t)] 



(3) 



For a discrete field, C is also equal to H[e{l~)]^ the Shan- 
non entropy of the local causal state^. C is the amount 
of information required to describe the behavior at that 
point, and equals the log of the effective number of causal 
states, i.e., of different distributions for the future. Com- 
plexity lies between disorder and order j4l|, |4j, 14^ , and 
(7 = both when the field is completely disordered (all 
values of x are independent) and completely ordered {x 
is constant). C grows when the field's dynamics be- 
come more flexible and intricate, and more information 
is needed to describe the behavior. 

The local complexity field C(r, t), is simply 
— logPr(s(r, t)) where s{f^t) is the local causal state 
at space-time point (r, t). The local complexity is the 
number of bits which would be required to encode the 
causal state at r, if we used the optimal (Shannon) 
coding scheme. Equivalently, it is the number of bits of 



information about /~(r, t) which are used to determine 
the causal state. 

It is appropriate, at this point, to take a step back 
and consider what we are doing. Why should we use the 
light-cone construction, as opposed to any other kind of 
localized predictor? Indeed, why use localized statistics 
at all, rather than global methods? Let us answer these 
in reverse order. The use of local predictors is partly a 
matter of interest — in studying coherent structures, we 
care essentially about spatial organization, and so global 
approaches, which would treat the system's sequence of 
configurations as one giant time series, simply don't tell 
us what we want to know. In part, too, the local ap- 
proach makes a virtue of necessity, because global pre- 
diction quickly becomes impractical for systems of any 
real size. The number of modes required by methods at- 
tempting global prediction, like Karhunen-Loeve decom- 
position, grows extensively with system volume jiM l4^. 
Global methods thus provide no advantage in terms of 
compression or accuracy. 

The use of light-cones for the local predictors, first sug- 
gested by 39],^ rather than some other shape, is moti- 
vated partly by physical considerations, and partly the 
nice formal features which follow from the shape, of which 
we will mention three [3^ . 

1. The light-cone causal states, while local statistics, 
do not lose any global predictive power. To be pre- 
cise, if we specify the causal state at each point in 
a spatial region, that array of states is itself a suf- 
ficient statistic for the future configuration of the 
region, even if the region is the entire lattice. 

2. The light-cone states can be found by a recursive 
filter. To illustrate what this means, consider two 
space-time points, (r, t) and (q^u) u > t. The 
state at each point is determined by the configu- 
ration in its past light-cone: s{f^t) = e{l~{f^t)), 
s{q^u) = e{l~{q^u)). The recursive-filtration prop- 
erty means that we can construct a function which 
will give us s{q^ u) as a function of s(r, t), plus the 
part of the past light-cone of (g, u) that is not vis- 
ible from (r , t) . Not only does this greatly simplify 
state estimation, it opens up powerful connections 
to the theory of two-dimensional automata 32]. 

3. The local causal states form a Markov random field, 
once again allowing very powerful analytical tech- 
niques to be employed which would not otherwise 
be available [50^ . 

In general, if we used some other shape than the light- 
cones, we would not retain any of these properties. 



2 The proof is as follows. /[e(/-);r] = if[e(r )] -if [e(/-)|r], the 
amount by which the uncertainty in is reduced by know- 

ing But for any discrete-valued function /, H[f{x)\x\ = 

0, because a function is certain, given its argument. Hence 
I[e{l-);l-] = H[e{l-)]. 



^ It seems that the use of light cones to analyze spatial stochastic 
processes was first suggested by Kolmogorov, who called them 
"causal sets", in the context of a model of crystallization |48| — 
see |49,|. 
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U ^ list of all pasts in random order 
Move the first past in U to a new state 
for each past in U 
noMatch ^ TRUE 

state ^ first state on the list of states 
while (noMatch and more states to check) 

noMatch ^ (Significant difference between 
past and state?) 

if (noMatch) 

state ^ next state on the list 
else 

Move past from U to state 
noMatch ^ FALSE 
if (noMatch) 

make a new state and move past into it from U 

FIG. 7: Algorithm for grouping past light-cones into esti- 
mated states 



B. A Reconstruction Algorithm 

We now sketch an algorithm to recover the causal 
states from data, and so estimate C. (See Fig. [7| for 
pseudo-code, Ref. js^ for details, and cf. 39].) At each 
time t, list the observed past and future light-cone con- 
figurations, and put the observed past configurations in 
some arbitrary order, {^^~}. (In practice, we must limit 
how far light-cones extend into the past or future.) For 
each past configuration estimate P^(/+|/~). We want 
to estimate the states, which ideally are groups of past 
cones with the same conditional distribution over future 
cone configurations. Not knowing the conditional distri- 
butions a priori, we must estimate them from data, and 
with finitely many samples, such estimates always have 
some error. Thus, we approximate the true causal states 
by clusters of past light-cones with similar distributions 
over future light-cones; the conditional distribution for a 
cluster is the weighted mean of those of its constituent 
past cones. Start by assigning the first past, to the 
first cluster. Thereafter, for each /~, go down the list of 
existing clusters and check whether P^(/+|/~) differs sig- 
nificantly from each cluster's distribution, as determined 
by a fixed-size ^^st. (We used a = 0.05 in all our 
calculations.) If the discrepancy is insignificant, add 
to the first matching cluster, updating the latter's dis- 
tribution. Make a new cluster if /~ does not match any 
existing cluster. Continue until every is assigned to 
some cluster. The clusters are then the estimated causal 
states at time t. Finally, obtain the probabilities of the 
different causal states from the empirical probabilities of 
their constituent past configurations, and calculate C(t). 
This procedure converges on the correct causal states as 
it gets more data, independent of the order of presenta- 
tion of the past light-cones, the ordering of the clusters, 
or the size a of the significance test 113 • Fo^" finite data, 
the order of presentation matters, but we finesse this by 
randomizing the order. 



C. Results on Elementary Cellular Automata 

In MBI we demonstrated the ability of local sensitivity 
filtering to find coherent structures in elementary cellu- 
lar automata. Here we apply local statistical complexity 
filtering to the same cellular automata, and compare the 
results of the two procedures. 

Figures Ofc and ^ show the local complexity fields 
C(r, t) of two of the classic rules, 110 and 54, respectively. 
Particles stand out clearly and distinctly on clean back- 
grounds. Note that we achieve this result by applying 
the same filter to both systems, even though their back- 
ground domains and particles are completely different. 
Were one to use, say, the conventional regular-language 
filter constructed in |^ for rule 54 on rule 110, it would 
produce nonsense. (Ref. |25^ discusses the domains, par- 
ticles and conventional filters of both these rules.) Note 
that the filter in 31] is hand-crafted, based on a detailed 
understanding of the CA's dynamics, whereas our filter, 
as we have said, is completely automatic and requires no 
human intervention. One might expect this generality 
to be paid for in a loss of resolving power, or missing 
system-specific features, but this does not appear to be 
the case^. For instance, the regular-language analysis of 
rule 54 31] identifies a subtle kind of particle, which con- 
sists of a phase shift in the spatially-periodic background 
domain. These phase shifts are hard to identify by eye, 
but show up very cleanly as particles in Figure^. 

For completeness, and further comparison with the lo- 
cal sensitivity, we include rules 22 and 146 filtered with 
respect with statistical complexity in Figures [Hfc andEJ:. 
Statistical complexity is also able to distinguish between 
the two CA configuration fields, even though the differ- 
ence is not readily apparent to the naked eye. 

Experimentally, we find that the depth of the past light 
cone is more relevant to proper filtering than the depth 
of the future light cone. (This may be related to the 
recursive estimation property of the causal states.) On 
elementary cellular automata, depth 2 is often sufficient, 
in the sense that further extensions of the cones do not 
change the states identified, but the results presented 
here use depth 3 for both future and past light cones. 
(Rule 41, not shown for reasons of space, required a past- 
depth of 5 in order to reach convergence.) 

While both our techniques make it easy to identify ob- 
jects like particles, even when they were previously hard 
to detect, they are very different filters, not just in their 
definition but also in their results. This can immediately 
be seen from the correlation coefficients^ of the sensitiv- 



This statement must be qualified by a recognition that we must 
supply the filter with enough data for it to find the right states. 
The learning rate of the state reconstruction algorithm is an 
important topic, beyond the scope of this paper. 

^ The correlation coefficient of statistics, pxv = — — — — is a 

dimensionless counterpart to the correlation function of statisti- 
cal physics, CxY — {-^^} ~ normalized so that its value 



ity and complexity fields — for rule 110, for instance, the 
correlation is a negligible 0.014. 

More abstractly, statistical complexity is a local quan- 
tity that is calculated using a global object, namely the 
probability distribution over causal states. Its accurate 
estimation thus needs a quite large number of cells, since, 
lacking an analytical form for that distribution, the latter 
must itself be estimated. A further consequence is that 
the separation between particles and other structures and 
the background tends to become cleaner over time, as the 
domains grow and the density of particles decays, sup- 
pressing the complexity of the former and raising that 
of the later. (This may be clearly seen in Figure 
Identification of structures from the complexity field is 
thus best undertaken after a (hopefully short) transient 
regime has passed, during which the local sensitivity fil- 
ter may be more useful. It is not, of course, necessary for 
the system to have reached a stationary regime in order 
to use the local complexity filter. 



IV. RESULTS ON CYCLIC CELLULAR 
AUTOMATA 

In this section we compare and contrast local sensi- 
tivity and statistical complexity with a more traditional 
order parameter approach in a spiral-forming CA model 
of an excitable medium. By applying all three analyses to 
the same 2+ ID CA, we will clarify the specific structural 
details emphasized by each method, and bring out the 
advantages of automatic filtering over approaches where 
the details must be put in by hand. 

Our model system in this section is the cyclic cellular 
automata (CCA) on a square lattice, which have been 
studied in considerable detail in the literature on spa- 
tial stochastic processes |^ [H^- In the general case, 
there are k colors, numbered from to z^: — 1. A cell of 
color k changes its color only if at least T of its neigh- 
bors are of color k -\- 1 mod k,^ in which case it assumes 
that color. In this paper, we confine ourselves to the spe- 
cial case of range-one neighborhoods, n = and T = 2. 
The behavior of the system, started from uniform ran- 
dom initial conditions, is illustrated in Figure [3 In [3^ . 
we demonstrated the ability of statistical complexity to 
quantify the extent to which the CA self-organizes, and 
the effect of changing parameter values on the degree of 
self-organization. Here, however, we are more interested 
in the patterns formed than in whether significant pat- 
tern formation is taking place, so we deal only with the 
most strongly self-organizing variant. 
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FIG. 8: (Color online.) A typical configuration of the spiral- 
forming CCA cyclic cellular automaton in the asymptotic 
regime (t = 200), illustrating the presence of coherent struc- 
tures (rotating spiral waves) dominating the lattice. 



A. Order Parameter and Effective Free Energy 

While, as remarked, cyclic CA have been extensively 
studied as models of excitable media, they have not (to 
the best of our knowledge) hitherto been examined with 
the tools of the order parameter approach. To better 
highlight the distinctive characteristics of the automatic 
filtering methods we propose, we first identify an order 
parameter and effective free energy for the present ver- 
sion of CCA. Our free energy will be a functional of the 
configuration field configuration alone, and is most ap- 
propriate for the long-run stationary distribution, rather 
than the initial transient stages, when the system is far 
from statistical equilibrium. Accordingly, we suppress 
time as an argument to the configuration field in what 
follows. 

As noted, with these parameter settings, the CA con- 
figuration field forms rotating spiral waves, which grow 
to engulf the entire lattice, with disordered domain walls 
at the boundaries between competing spiral cores. The 
mechanisms by which spirals form and grow are fairly 
well understood, and topological arguments 52] pick out 
the key role of both conservation of winding number and 
of the spiral cores in this process. Accordingly we con- 
sider CCA as a kind of discretized XY model, similar to a 
clock model 4, §3.6.3]^, and, as in such models generally. 



lies between —1 and +1, and is zero when X and Y are linearly 
unrelated. 



^ A class of cellular automata very similar to CCA are treated as 
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the appropriate order parameter has two components |j . 
Here the crucial observation is that the ground state con- 
sists of plane waves, with stripes of cells of constant color 
extending perpendicular to the direction of propagation 
of the wave. The order parameter is the local normalized 
wave- vector. 

^(x, y) = ^^(x, y)x + ^y{x, y)y 

defined such that it takes a one of four values ^ = 
^(±1, ±1) in each of the 4 domains surrounding a spi- 
ral core. The exact definition of the wave vector in terms 
of the states of the CCA configuration field is given in 
appendix |B| This wave vector can be used to define a 
phase for the spiral wave at each lattice site. 



<l>(x, y) = tan 



-1 %{x^y) 

^x{x,y) 



The local free energy is then calculated as the discretized 
version of 



F{x,y) = iW^x,y)Y 



(4) 



We show the free energy of the CA configuration field 
in Figure El Note the increased free energy caused by 
topological defects such as the domain walls and spiral 
cores. 



B. Local Sensitivity and Statistical Complexity for 

CCA 



Filtering by sensitivity and complexity, which we 
demonstrated above for elementary cellular automata, 
works without modification on cyclic CA. As might be 
expected, the two filters reveal different, but compatible, 
aspects of the system. 
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FIG. 9: Free energy per site of the CCA configuration field 
given in Figure |H| Darker cells have higher energy values. 
Note the increased free energy at the topological defects, e.g. 
the spiral cores and domain walls. 



antiferromagnetic Potts models in |53| . We experimented with 
such an order parameter, but the results were poor. 



FIG. 10: Local sensitivity (p = 0, / = 5) for the cyclic cellular 
automaton, using the CCA configuration field of Figure |H1 

Figure irni shows the local sensitivity field for the same 
configuration as in Figure |H1 with a one-cell perturbation 
range {p = 0). Figure ITTl on a smaller field, shows that 
with p = 1 we get qualitatively similar results. The spi- 
ral cores are easily spotted in both figures, demonstrating 
that they are autonomous objects. There is no clear dif- 
ference between spiral angles (the horizontal and vertical 
lines) and spiral domains. Spiral boundaries are white 
(lowest sensitivity), because any perturbation here will 
be quickly erased under the pressure of the radiating spi- 
ral cores (cf. Calculating local sensitivity for 2+ ID 
automata is very slow: each cell has 8 neighbors and there 
are 4 states, which makes 4^ — 1 possible perturbations. 
While sampling on a random subset of all possible per- 
turbations would be faster, the resulting approximation 
error would have to be carefully determined. 

We show the CCA field filtered with respect to statis- 
tical complexity in Figure El Comparison with Figure 
El shows that causal states and order parameter analysis 
yields to the same information. The spiral cores are still 
among the most complex areas. Here spiral boundaries 
are also complex: under statistical complexity filtering, 
this is because some prediction is possible, but requires 
much information; under free energy analysis, this is due 
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FIG. 11: Local sensitivity (a) for the cyclic cellular automa- 
ton, calculated with p = 1, / = 1, and corresponding config- 
uration (b). 



to high phase differences. 



C. Comparison of Free Energy, Sensitivity and 
Complexity 

The three methods described here — the free energy, 
the local sensitivity and the statistical complexity — all 
uncover significant structures in the spatiotemporal dy- 
namics of the CA. However, the methods are not equiva- 
lent, and emphasize different aspects of those structures, 
as one might expect of different filters. 

Qualitatively, we can see the difference in the order in 
which different features are ranked. For sensitivity, the 
ordering is particles (spiral cores), domains, and then 
domain walls last. For both complexity and free energy, 
the order is by contrast domain walls, particles and do- 
mains. A quantitative calculation of the correlation coef- 
ficients between the different filtered fields confirms this 
qualitative impression. The most strongly correlated — 
as one might expect by comparing Figure El with Figure 
[n\ — are the statistical complexity and the free energy, 
p = 0.690; we discuss the roots of this strong correla- 
tion below. The complexity and the sensitivity, however, 
are almost completely uncorr elated, just as we observed 
in the case of rule 110 — the correlation coefficient calcu- 
lated from sample data, p = —0.006, is not materially 
different from 0. Unsurprisingly, there is also no impor- 
tant correlation between the sensitivity and the free en- 
ergy {p = —0.008). Despite the fact that sensitivity and 
complexity are nearly orthogonal, plotting complexity as 
a function of sensitivity (not shown) reveals an interest- 
ing relationship: as sensitivity increases, the minimum 
value of complexity rises, though the converse is not true 
(minimally sensitive points are found at all values of com- 
plexity). 

Sensitivity looks for regions that are unstable to per- 
turbation; because of this instability, prediction requires 
very precise discriminations among histories. Moreover, 
localized unstable regions are presumably rare, and the 
corresponding states uncommon. (Recall that the com- 
plexity is both H[S] and I[S;X~].) Thus the spiral cores 
are rare, sensitive and complex. On the other hand, un- 
common, complex structures need not be unstable. Do- 




FIG. 12: Local statistical complexity of the CCA configu- 
ration field in Figure |H1 calculated with depth 1 light-cones. 
Darker cells have higher values for the field. This figure should 
be compared to both Figures El and Note that free energy 
and statistical complexity are strongly correlated, while lo- 
cal sensitivity and statistical complexity are not because they 
emphasize different aspects of the coherent structure. 



main walls, for instance, are comparatively rare, and re- 
quire considerable historical information for their predic- 
tion, essentially because they are regions where the evolu- 
tion of the phase is hard to determine; this in turns means 
that their causal states need more bits to specify. They 
are, however, very stable, because they have consider- 
able spatial extent, and to destroy or move a wall implies 
changing the domains on either side. The domains them- 
selves, while only minimally complex, are more sensitive 
than the domain walls, because perturbations there can 
create localized defects. 

The relationship between free energy and complexity, 
while strong and at first sight clear, is actually a bit 
tricky. The easy, but flawed, argument for why effective 
free energy and complexity should be correlated runs as 
follows. The effective free energy relies on the assumption 
that the microscopic dynamics tend to create ordered re- 
gions, where "ordered" is understood in the sense given 
by the order parameter. All departures from the ordered 
state are unlikely, and large departures are exponentially 
unlikely. (This is not tautologically true, but when it is, 
we can define an effective free energy.) These departures 
take the form of topological defects, like spiral cores and 
domain walls, which pay an energetic penalty. The ex- 
act size of this penalty, and the resulting rarity of these 
features, will depend on the details of the microscopic dy- 
namics. But, precisely because they are rare, the causal 
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states corresponding to them should be rare too, and so 
have high local complexity. 

The flaw is that the causal state at a point is a func- 
tion of the configuration in its past light cone, while the 
order parameter and the free energy are functions of the 
immediate neighborhood configuration. Since the two fil- 
ters use different inputs, their strong correlation is not 
trivial. However, under fairly general conditions, the lo- 
cal causal states form a Markov random field, and the 
Gibbs-Markov theorem tells us that a random field is 
Markovian if and only if there is an additive effective 
free energy which gives the probability of configurations 
[S^. Conversely, it can be argued that a complete set of 
thermodynamic macro-variables should be equivalent to 
the causal state of the system, and should be Markovian 
in time, so that they can be calculated entirely from the 
present configuration Thus, there should, in fact, 

be a relationship between the complexity and the free 
energy, if we have identified the proper order parame- 
ter; the strong correlation we find between the two fields 
suggests that we have done so. The fact that the corre- 
lation, while strong, is not perfect could be either due to 
our definition of the order parameter being slightly off, 
or due to finite-sample errors in the identification of the 
causal states. 



V. CONCLUSIONS 

We have introduced two complementary filtering meth- 
ods for spatially extended dynamical systems. Local sen- 
sitivity calculates the degree to which local perturbations 
alter the system, and picks out autonomous features. Lo- 
cal statistical complexity calculates the amount of his- 
torical information required for optimal prediction and 
identifies the most highly organized features. We em- 
phasize that "organized" is not equivalent to ordered, 
because the latter corresponds to low entropy which is 
not the same as high complexity or organization |3^. A 
regular lattice is highly ordered, but not very organized. 
On the other hand, high complexity is not equivalent 
to high entropy. A random field has high entropy but 
low complexity. Complexity, as we said earlier, lies be- 
tween order and randomness, and is a refiection of the 
probabilistic properties of the system. In contrast, sensi- 
tivity is a measure of the system's dynamical properties. 
The two are linked through ergodic theory, but they re- 
main distinct. Both sensitivity and complexity pick out 
spatiotemporal coherent structures, and the structures 



^ Eric Smith (personal communication) has pointed out that this 
argument suggests a relationship between the time depth needed 
to identify the causal states, and the spatial range needed to cal- 
culate the free energy. However, pinning down that relationship 
would require careful treatment of many mathematical issues re- 
garding sofic shifts, long-range correlations in Markovian fields, 
etc., beyond the scope of this paper. 



they identify match those known from previous, more ad 
hoc approaches, whether based on regular languages (as 
in ECA) or order parameters and topological considera- 
tions (as in CCA). In no case, however, is prior knowledge 
about the system or its coherent structures used in con- 
structing our filters; at most we have tuned calculational 
parameters, in a way akin to adjusting a microscope until 
the image comes into focus. Ideally, one would use both 
sensitivity and complexity filtering, because they provide 
distinct kinds of information about the system, but we 
suspect complexity will be much easier to calculate from 
empirical data since it only requires observation and not 
perturbation. 

Many theoretical and mathematical questions present 
themselves about both filtering methods. We have men- 
tioned some in passing; here we wish to highlight just a 
few. On the local sensitivity method: What is the exact 
relationship between the perturbation range p and the 
spatial scale of identified structures? How much error 
would be introduced by averaging over a random sub- 
set of perturbations, rather than an exhaustive enumera- 
tion? Can one identify a typical lifespan for a perturba- 
tion, after which it is erased by its surroundings, and if 
so is this lifespan related to /? (This last is presumably 
related to dynamical mixing properties.) On the local 
complexity method: what quantitative factors relate the 
volume of data available to the error in our estimates of 
the causal states and so of the complexity? Can we use 
causal states to give algebraic/automata-theoretic defi- 
nitions of "domain" and "particle", like those in the ID 
case, without entangling ourselves in the difficulties of 
higher-dimensional languages? (Cf. 36, §10.5].) When 
we do have partial knowledge of the correct pattern ba- 
sis, who can we use this to hasten the identification of the 
local causal states? Perhaps most ambitiously, could one 
reverse engineer a good order parameter from the local 
causal state field and its transition structure? 

As very large, high-dimensional data sets become in- 
creasingly common and important in science, human per- 
ception will become increasingly inadequate to the task 
of identifying appropriate patterns [s^, IE3 • desir- 
able, therefore, to move towards more automatic filter- 
ing methods, and automatic ways of detecting coherent 
objects. Because our filtering methods do not presup- 
pose any prior knowledge or require human insight about 
what the right structures are, they should work generi- 
cally, across systems of highly varying nature and dy- 
namics. This is a hypothesis rather than a theorem, but 
it can be tested simply by applying our filtering meth- 
ods to a wide range of systems with known coherent 
structures — and to others where the appropriate struc- 
tures are not known. In the latter cases, the test of our 
methods will be whether the structures they identify can 
be used to frame interesting and insightful higher-level 
models about the dynamics and functions of the systems 
involved (cf. IHil). 
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APPENDIX A: ATTEMPTS TO ADAPT 
LYAPUNOV EXPONENTS TO CELLULAR 
AUTOMATA 

There have been at least three attempts to define Lya- 
punov exponents for cellular automata. We review them 
in order of priority, concluding that none of them is alto- 
gether suited to the aims of this paper. All are based on 
measuring the rate at which small perturbations cause a 
divergence from the original trajectory. 

Shereshevsky 59] defines a Lyapunov exponent for 
one-dimensional cellular automata, as the velocity of 
propagation of the edge of the difference plume. More 
exactly, he defines two Lyapunov exponents, for the ve- 
locities of the left and right edges. By assuming an in- 
variant measure over configurations exists and has certain 
properties, he is able to show that these velocities have 
reasonable long-time limits. (Tisseur shows the limit 
exists with somewhat weaker j conditions on the invariant 
measure.) To show that the limit is uniform over the lat- 
tice, he invokes a further spatial ergodicity property. 

We, of course, would like things not to be uniform, so 
the fact that we don't have measures which are so nicely 
ergodic and stationary should not trouble us. But this 
definition allows a cell which only makes a difference to 
one other cell in the future to count as highly infiuential, 
provided that said cell moves very rapidly. (Consider the 
shift rule 170 "copy your neighbor to your left" .) Looking 
at the area of the difference plume (i.e., the number of 
differing cells between the original configuration and the 
perturbed configuration) seems far more reasonable. 

Bagnoli et al. 61] come closest to the way we define 
local sensitivity in section^ by looking at the total num- 
ber of "defects" , here meaning the number of cells which 
are different between the original and the perturbed con- 
figuration. They proceed as follows. They perturb a sin- 
gle site of the lattice, and iterate forward one time step, 
so that the the difference plume now embraces m sites. 
They then create m copies of the lattice, each identical 
to the time-evolved unperturbed configuration, except at 
one of the m sites in the difference plume. They then re- 
peat this procedure, accumulating more and more copies 
of the lattice, each of which differs from the unperturbed 
trajectory only at a single site. Suppose the initial per- 
turbation was applied at cell xq and time to. Then for 
each site x and time t > to, the number of copies of the 
system which have a defect at x at time t is N^^^toix^t)^ 



and N^o,to{^) = J2x^i^^^)' ^^o,to(^) can grow expo- 
nentially, and Bagnoli et al. define the local Lyapunov 
exponent to be that exponential growth rate. 

This is far from anything that could be called a Lya- 
punov exponent in the strict sense of the term. Fur- 
thermore, to avoid the extremely involved (exponential) 
calculation which their definition implies, they make use 
of a kind of derivative, as in the continuous-system defini- 
tion of the Lyapunov exponent. But these derivative are 
only defined for Boolean- valued rules, which makes them 
useless for, e.g., cyclic cellular automata. The direct cal- 
culation could be somewhat simplified for deterministic 
cellular automata by means of dynamic programming, 
but it would still be a hard calculation to obtain a num- 
ber whose physical significance is unclear. 

Finally, Urias et al. 62] define a quantity which is 
a Lyapunov exponent (and doesn't require binary val- 
ues), but only under a very specific metric, with no clear 
physical interpretation; it is not clear how, if at all, their 
metric could be extended to higher-dimensional cellular 
automata. The ultimate result they get is that the Lya- 
punov exponent is just the maximum velocity at which 
the envelope of the difference region spreads — taking the 
maximum over all possible semi-infinite (not single-point) 
perturbations. 

In this paper we wished not to measure a single Lya- 
punov exponent for the entire system, but rather the local 
(in both space and time) effects of perturbations. This 
allowed us to determine the degree to which different lo- 
cal structures of the CA are autonomous. 



APPENDIX B: ORDER PARAMETER AND 
FREE ENERGY FOR SPIRAL-FORMING 
CYCLIC CELLULAR AUTOMATA 

A classic method of describing the equilibrium ordered 
phases of a system is to define an appropriate order pa- 
rameter for each phase and an associated free energy 
determined from symmetry considerations. In this ap- 
pendix we will perform this analysis for the spiral cellular 
automaton in its equilibrium state. By "equilibrium" we 
mean the cellular automaton's long time behavior after 
the domain walls and spiral cores have stabilized. The 
transient period of the cellular automaton corresponds 
to the phase transition which is not described by the 
order parameter formalism. Examination of the cellu- 
lar automaton field (e.g.. Figure [5|) reveals that, while it 
certain possesses order, it is also highly frustrated and 
exhibits many topological defects, most notably the vor- 
tices at the spiral centers and the domain walls between 
adjacent spirals. However the order parameter is defined, 
we expect the free energy to increase at these defects. 

We begin by defining an appropriate order parameter 
on the discrete lattice. We take our cue from the fact that 
for a continuous field, the spirals would consist of con- 
centric rings possessing two dimensional rotational sym- 
metry. We therefore construct a discretized XY model 
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order parameter and free energy. In the XY model the 
Ginsburg-Landau free energy density at a given lattice 
site is: 

F{x,y) = a{V^x,y)f (Bl) 

a is a positive constant which we arbitrarily set to one 
and henceforth ignore. <l> is a localized phase defined at 
each lattice site. In practice we calculate the Laplacian 
using each lattice site's nearest neighbors. 

Fix,y) = ^{i^ix + l,y)-^ix,y)f 

+ {^{x,y)-'^{x-l,y)f 
+ (*(a;,y + l)-*(a;,y))2 
+ {^{x,y)-^{x,y-l)f} (B2) 

$ is defined in terms of the states of the configuration 
field a e {0,1,2,3}. Noting that the spiral cellular 
automaton field consists mainly of plane waves cycling 
through the four colors and traveling mainly along the 
diagonals, although at times along the x and y axes, it 
seems reasonable to define the spiral cellular automaton 
phase ^ as a function of the direction in which the wave 
is traveling. We define a normalized two component wave 
vector at each lattice site (x, y) 

^{x,y) = ^ {^^{x,y)x^^y{x,y)y) 

= i^x.^y) (B3) 

where and are the components of the normalized 
vector and the unnormalized components and are 
defined as follows. 

^^xix.y) = ^ g{a{x^i,y),a{x,y)) (B4) 

+ ^ 9i^ix + hy^j)^cr{x,y)) 

i=±ij=±i 



"^yi^^y) = "Y 9{(^{x,y ^i),cr{x,y)) (B5) 
i=±i 

+^ Yl 9[(^{^^hy^ j)^(^{x^y)) 
i=±i,j=±i 

Note that the first sum is over nearest neighbors in the 
X (resp. y) direction and the second sum is over all next 
nearest neighbors. The factor 1/2 comes from the cor- 
responding 1/2 in the definition of the order parameter. 
The function g is defined as 

g{(7l,(72) = +1 if CTi = (72 + 1 
= — 1 if (Ti = <J2 — 1 

= otherwise (B6) 

It is to be understood that the above definition of g in- 
volves addition and subtraction modulo 4. 

The phase <l> at a given lattice site is defined simply as 

$ = tan-^ ^ (B7) 

Thus we see that in the upper right quadrant of a spiral, 
^ = (^, ^) and <l> = ^, whereas in the upper left quad- 
rant ^ (^, and = Thus there is a phase 
gradient along the boundary between these two quad- 
rants of the spiral and an increase of free energy here. 
This can be seen in Figure El (The localization of the 
free energy increase to the boundary between quadrants 
is, of course, an artifact of the discretized spatial lattice. 
In the continuum limit, the gradient energy would be 
distributed evenly around the central vortex.) The other 
locations where one expects to see increased free energy 
are in the vortex cores, and at the boundaries between 
spirals. Our free energy captures both of these effects 
(Figure 01). 
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